restore,'grid.sav'
if !d.name eq 'PS' then begin
    device,xsize=25,ysize=6,yoffset=3
    !p.thick=1 & !x.thick=1 & !y.thick=1
end
@rhoprof

restore,'velocity.sav'
restore,'alpha_eta.sav'
alpha=transpose(khel)
nr=80L
nth=80L
tth=!pi/2.+congrid(th[0:31],nth,cubic=-0.5)*!pi/(180.)
rr=congrid(r,nr,cubic=-0.5)
x=rr#sin(tth)
y=rr#cos(tth)
oo=congrid(omega[0:40,4:31],nr,nth,cubic=-0.5)
ur=congrid(ww[0:44,4:31],nr,nth,cubic=-0.5)
uth=congrid(vv[0:44,4:31],nr,nth,cubic=-0.5)
aa=congrid(alpha[0:44,4:31],nr,nth,cubic=-0.5)
rhoxy=congrid(rh2d,nr,nth,cubic=-0.5)
;;;;;;;;;;;;;;;;;;;;;;;;;;;
;normalization coefficients
oo0=2.*!pi/(28.*24.*3600.)
eta0=1.e12
aa0=10.
	
;;;;;;;;;;;;;;;;;;;;;;;;;;
etat=smooth(etat,2)
eta=congrid(etat,nr)/eta0
aa=aa/aa0

oo=2.*!pi*smooth(oo,3)*1.e-9/oo0
ur=smooth(ur,3)
uth=smooth(uth,3)
aa=smooth(aa,3)
;eta=smooth(eta,3)
umax=max([abs(uth),abs(ur)])
uth=uth/umax
ur=ur/umax

doodr=fltarr(nr,nth)
durdr=fltarr(nr,nth)
duthdth=fltarr(nr,nth)
;r derivatives
for j=0,nth-1 do begin
 doodr[*,j]=xder_curl(oo[*,j],rr)
; doodr[*,j]=deriv(rr,oo[*,j])
 durdr[*,j]=xder_curl(ur[*,j],rr)
endfor
detadr=smooth(xder_curl(eta,rr),2)
;detadr=smooth(deriv(rr,eta),2)

doodth=fltarr(nr,nth)
for i=0,nr-1 do begin
 doodth[i,*]=xder_curl(reform(oo[i,*]),tth)
 duthdth[i,*]=xder_curl(reform(uth[i,*]),tth)
endfor

fo='(F9.6,1x,F9.6,1x,9(F11.6,1x))'  ; ,2x,8(E11.4,2X))'
openw,1,'dynamo.dat'
for i=0,nr-1 do begin
 for j=0,nth-1 do begin
  print,rr[i],tth[j],doodth[i,j],doodr[i,j],aa[i,j],ur[i,j],durdr[i,j],uth[i,j],duthdth[i,j],eta[i],detadr[i],fo=fo
  printf,1,rr[i],tth[j],doodth[i,j],doodr[i,j],aa[i,j],ur[i,j],durdr[i,j],uth[i,j],duthdth[i,j],eta[i],detadr[i],fo=fo
 endfor
endfor
; dynamo numbers
dr=r[1]-r[0]
maxshear=max(abs(doodr))
print,'DYNAMO NUMBERS '
rrr=6.96e10
cw=oo0*(rrr)^2/max(etat)
ca=100.*max(abs(alpha))*rrr/max(etat)
cu=100.*umax*rrr/max(etat)
print,'cw = ',cw
print,'ca = ',ca
print,'cu = ',cu
print,'Re= ',max(etat)/min(etat)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; plotting
!p.multi=[0,4,1] & !p.charsize=1.7 
rsinth=fltarr(nr,nth)
r2sinth=fltarr(nr,nth)
for k=0, nr-1 do begin
  for j=0,nth-1 do begin
    rsinth[k,j]=rrr*rr[k]*sin(tth[j])
    r2sinth[k,j]=(rrr*rr[k])^2*sin(tth[j])
  endfor
endfor

rhov=uth*rhoxy*rsinth
rhow=ur*rhoxy*r2sinth
stream1=fltarr(nr,nth)
stream2=fltarr(nr,nth)
for j=0,nth-1 do begin
  for k=0,nr-3 do begin
    kk=indgen(k+2)
    stream1[k,j] = int_tabulated(rr[kk],rhov[kk,j],/double)
  endfor
endfor
  
for k=0,nr-1 do begin
  for j=0,nth-3 do begin
    jj=indgen(j+2)
    stream2[k,j] = int_tabulated(tth[jj],rhow[k,jj],/double)
  endfor
endfor
thg=where(tth gt 0)
stream=smooth(stream1+stream2,2)

nls=10
lev_stream = 2.*max(abs(stream[*,thg]/rsinth[*,thg]))*indgen(nls)/float(nls-1)$
             -max(abs(stream[*,thg]/rsinth[*,thg]))
contour,oo,x,y,/fi,nl=64,/iso,title='!7X!6'
contour,stream[*,thg]/rsinth[*,thg],x[*,thg],y[*,thg],/over,lev=lev_stream, $
        max_value=0.,min_value=min(lev_stream),c_thick=2,c_linestyle=2,color=255
contour,stream[*,thg]/rsinth[*,thg],x[*,thg],y[*,thg],/over,lev=lev_stream, $
        max_value=max(lev_stream),min_value=0.,c_linestyle=0,c_thick=3,color=255
contour,doodr,x,y,/fi,nl=64,/iso,title='!4d!7X!6/!4d!8r!6'
contour,doodth,x,y,/fi,nl=64,/iso,title='!4d!7X!6/!4d!4h!6'
contour,aa,x,y,/fi,nl=64,/iso,title='!4a!6'
!p.multi=0
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
close,1
print,'umax= ',umax
print,'UTHmax = ',minmax(uth)
print,'Ur = ',minmax(ur)
print,'minmax alpha = ',minmax(aa)
print,'minmax eta = ',minmax(etat)


end
